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Error  Analysis  of  Algorithms  for  Evaluating 
Bernstein- Bezier-Type  Multivariate  Polynomials 

J.  M.  Pena 


Abstract.  In  Computer  Aided  Geometric  Design,  the  Bernstein-Bezier 
form  is  the  usual  way  to  store  a polynomial  defined  on  a triangle.  We  per- 
form backward  and  forward  error  analysis  of  the  de  Casteljau  algorithm 
and  of  the  algorithm  proposed  by  Schumaker  and  Volk  for  evaluating  such 
polynomials.  The  obtained  results  are  also  compared  with  the  correspond- 
ing results  for  the  bivariate  Horner  algorithm. 


§1.  Introduction 

In  Computer  Aided  Geometric  Design,  multivariate  polynomials  defined  on  a 
triangle  are  usually  stored  in  the  Bernstein-Bezier  form,  and  can  be  evaluated 
by  the  de  Casteljau  algorithm.  In  [8]  a modified  Bernstein-Bezier  represents 
tion  of  polynomials  was  introduced,  along  with  an  algorithm  for  its  evaluation, 
which  was  called  the  VS  algorithm.  This  algorithm  for  the  evaluation  of  mul- 
tivariate polynomials  is  expressed  in  terms  of  nested  multiplications,  and  is 
more  efficient  than  the  de  Casteljau  algorithm. 

Error  analysis  of  the  de  Casteljau  algorithm  for  univariate  polynomials 
was  considered  in  [2]  and  [4] . This  paper  is  devoted  to  backward  and  forward 
error  analysis  of  the  de  Casteljau  and  VS  algorithms  for  bivariate  polynomi- 
als. On  the  other  hand,  the  error  analysis  of  the  Horner  algorithm  for  the 
evaluation  of  univariate  polynomials  has  been  extensively  studied  in  the  liter- 
ature. In  fact,  backward  and  forward  error  analysis  of  (univariate)  Horner’s 
rule  was  already  performed  by  Wilkinson  in  [11],  pp.  36-37  and  49-50.  Other 
approaches  to  this  problem  can  be  found  in  [5,9,10,12]  (see  more  references  in 
[3]).  An  error  analysis  of  the  multivariate  Horner  algorithm  has  been  given  in 
[7].  We  also  compare  our  results  with  those  corresponding  to  the  multivariate 
Horner  algorithm. 

The  paper  is  organized  as  follows.  Section  2 introduces  basic  concepts, 
notation,  and  auxiliary  results.  In  Section  3 we  carry  out  the  mentioned 
error  analysis  of  the  algorithms.  Finally,  we  summarize  in  Section  4 the  main 
conclusions  and  the  advantages  of  VS  algorithm  in  this  context,  taking  into 
account  computational  cost  and  forward  error  analysis. 
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§2.  Basic  Notations  and  Auxiliary  Results 

Let  us  now  introduce  some  standard  notation  in  error  analysis.  Given  aeR, 
the  computed  element  in  floating  point  arithmetic  will  be  denoted  by  either 
fl(a)  or  by  a.  As  usual,  to  investigate  the  effect  of  rounding  errors  when 
working  with  floating  point  arithmetic,  we  use  the  model 


fl(aopi)  = (aop6)(l  + 6),  |<5|  < u,  (1) 

with  u the  unit  roundoff  and  op  any  of  the  elementary  operations  +,  — , x,  / 
(see  [3,  p.  44]  for  more  details).  Given  k S N0  such  that  ku  < 1,  let  us  define 


7 k ■= 


ku 

1 — ku 


(2) 


In  our  error  analysis  we  shall  deal  with  quantities  such  that  their  absolute 
value  is  bounded  above  by  7*.  Following  [3],  we  denote  such  quantities  by 
6k  and  take  into  account  that  by  Lemmas  3.3  and  3.1  of  [3],  the  following 
properties  hold: 

(1  + 0*)(1  + Oj)  = l + 6k+j,  (3) 

and  if  p,  = ±1,  |5<|  < u (i  = 1, . . . , k)  then 


k 

n(i+*r=i+fc. 

*= 1 


(4) 


In  considering  the  computed  solution  of  a problem,  one  can  try  to  find  the 
data  for  which  this  computed  solution  is  the  exact  solution.  Backward  error 
analysis  measures  how  far  these  data  are  from  the  original  data  of  the  prob- 
lem. So,  backward  error  analysis  interprets  rounding  errors  as  perturbations 
in  the  data.  In  contrast,  forward  error  analysis  measures  how  far  the  computed 
solution  is  from  the  exact  solution.  Therefore,  in  our  evaluation  problem,  if 
f(x)  = is  the  computed  evaluation  (instead  of  the  exact  evalua- 

tion f(x)  = X]I*=o  °iU >(x))>  we  say  that  the  relative  backward  error  is  bounded 
above  by  e if 


Then  we  can  bound  the  forward  error  by 


\f(x)-f{x)\  <eJ2\ciUi{x)\. 

i= 0 


The  number 

n 

Cu(f{x))  ■■=  '52  |ciUi(x)|,  (5) 

i=0 

measures  the  stability  in  the  evaluation  of  a function  with  respect  to  perturba- 
tions of  the  coefficients,  and  is  called  the  condition  number  for  the  evaluation  of 
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/( x)  with  the  basis  u (see  [1,2, 4, 6, 7]).  Let  us  observe  that  Cu(f(x))  depends 
on  the  basis  u,  on  the  function  /,  and  on  the  point  x.  If  we  assume  that  the 
basis  is  formed  by  nonnegative  functions,  (5)  can  be  written  as 


Cu(f(x))  = \ci\ui{x). 

i=0 


In  conclusion,  we  can  bound  the  forward  error  by 

\f(x)-f(x)\  <eCu(f(x)), 


(6) 

(7) 


which  is  a particular  case  of  the  classical  formula  (Forward  error)  < (Backward 
error)  x (Condition  number). 

If  we  assume  that  f(x)  ^ 0 we  can  also  define  the  relative  condition 
number  by 


cu(f(x)) 


Cu(f(x)) 

\m\  ' 


The  following  relative  forward  error  bound,  analogous  to  (7),  can  be  derived: 


l/» -/(*)! 
!/(*)! 


< ecu(f(x)). 


The  following  result  (which  was  obtained  for  polynomials  in  [1])  allows 
us  to  compare  the  condition  number  of  two  bases  of  nonnegative  functions  in 
a space  when  they  are  related  by  a nonnegative  matrix: 


Lemma  1.  Let  It  be  a finite  dimensional  vector  space  of  functions  defined  on 
SI  C 1RS.  Let  u = (mq,  , un),  v = (no, ... , vn)  be  two  bases  of  nonnegative 
functions  ofU  such  that 


(vo,  • • • i ^n)  — (^0j  • • • j 

where  A = (aij)o<i,j<n  is  a nonnegative  matrix.  Then  Cu(f(x))  < Cv(f(x)) 
for  each  function  f € U evaluated  at  every  x £ SI. 

Proof:  Given  / € It,  it  can  be  written  as 

n n / n \ 

f(x)  = YjWiW  = «<(*)•  (8) 

9=0  i=0  \9=0  / 

Then,  by  (8)  and  the  nonnegativity  of  u,  v and  A,  we  deduce  that 

C«(/(^))  = ElEMi9)|k(t)| 

i=0  9=0 

n n n / n \ 

— |c9|aj9Uj(a:)  = |cg|  I aiqUi(x)  I (9) 

1=0  9=0  9=0  \l=0  / 

n 

= J2  w Mx)i = <?„(/(*)) 

9=0 

for  each  function  / eU  evaluated  at  every  x € □ 
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§3.  Backward  and  Forward  Error  Analysis  of  the  Algorithms 

Given  a triangle  T with  vertices  P,  Q,  R in  the  plane  and  a point  X £ T,  let 
( r , s,t)  be  the  barycentric  coordinates  of  X: 

X = rP  + sQ  + tR,  r + s + t = 1. 

Let  nrf(T)  space  of  polynomials  of  total  degree  d defined  on  T.  Then  any 
polynomial  p £ II d(T)  can  be  written  as 

d i 

p{r,  s,  t)  = Y.  S ’ Oi  (10) 

t=0  j= 0 

where 

B?tJtk(r,8,t)  = i + j + k = d (11) 

are  the  Bernstein  polynomials.  Then  (10)  is  called  the  Bernstein-Bezier  repre- 
sentation of  p. 

Let  us  now  recall  the  algorithm  of  de  Casteljau  to  evaluate  the  polynomial 
p at  a given  ( r,s,t ).  We  denote  by  :=  bitjtk  for  all  i + j + k = d. 


The  previous  algorithm  leads  to  the  evaluation  p(r,s,t ) = ci  o ■ ^ re“ 
quires  d(d  + l)(d  + 2)/2  multiplications. 

The  following  result  provides  backward  error  analysis  of  the  de  Casteljau 
algorithm  for  the  evaluation  of  bivariate  polynomials. 

Theorem  2.  Let  us  consider  the  algorithm  of  de  Casteljau  in  (10)  and  let 
us  assume  that  Mu  < 1.  Then  the  computed  value  p(r,  s,t)  = fl(p(r,  s,t)) 
satisfies 

d i 

P(r,s,t)  =Y^Y^bd-Li-idBd-i,i-jAr’s't')’  (12) 

i=0  j= 0 


where 


I bi,j,k  bi,j,k  I ^ 
— ii — i — - 73d 


(13) 
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hk 

ud—i—kti—j,j 

- [fi(rbd-i-k+l,i-j,j)  + ft(sbd-i-k,i-j+l,j  + tbd-i-k,i-j,j+d]  (*  + *>) 

= [(rfed-i-fc+i,j_j,j)(l  + <$i)  + (&(sbkdz]_kj_j+lj) 

+ + ^2)]  (1  + 60) 

= [(rfod-i-fc+1,i_jj)(  1 + <5i)  + ((s^Ii-fc,i-j+1,j)(l  + ^3) 

+ (^d-»1-fc,i-i,j+i)(1  + ^))(1  + ^2)]  (1  + 60), 


where  |5j|,  i = 0, . . . , 4,  are  numbers  less  than  or  equal  to  the  unit  roundoff  u. 
Then  by  (4)  we  can  write 


bd-i-k,i-j,j  - (ri>d-i-k+l,i-j,j)0-  + 02 ) + (s&d-t1-fc,i-i+l,j)(1  + 0s) 

+ (tbkdz]_k>i_jJ+1)(l  + 03). 

Iterating  the  previous  argument  for  k = d,  d — 1, . . . , 1 and  comparing  the 
final  expression  with  (10),  we  deduce  (12)  and  (13).  □ 

By  applying  the  previous  result  together  with  formula  (7),  which  relates 
backward  and  forward  error,  we  can  derive  the  following  corollary  on  the 
forward  error  of  the  de  Casteljau  algorithm. 

Corollary  3.  Suppose  the  hypotheses  of  Theorem  2 hold.  Then  the  com- 
puted value  p(r,  s,  t)  = fl(p(r,  s,  t))  given  by  the  de  Casteljau  algorithm  satis- 
fies 

\p(r,s,t)  -p(r:s,t)\  < yidCB{p{r,s,tj),  (14) 

where  B is  the  Bernstein  basis. 

A modified  Bernstein-Bezier  representation  of  the  polynomial  p of  (10) 
was  considered  in  [8], 


P( 


r,  S,  t)  — ^ ^ ) Cd—i,i—j,j 

i=0  j= 0 


„d—i  A 
T S 


(15) 


where  the  new  coefficients  are  related  with  those  of  (10)  by 


Cd—  i. 


d\ 


(d  — i)\(i  — j)\j\ 


■ibd—i,i—j,j:  0<j,i<d. 


The  following  algorithm  (which  will  be  called  the  VS  algorithm)  was  in- 
troduced in  [8]  by  Schumaker  and  Volk  to  evaluate  a polynomial  p in  the 
modified  Bernstein-Bezier  form  (15).  In  contrast  with  de  Casteljau  algorithm, 
this  algorithm  is  expressed  in  terms  of  nested  multiplications.  This  version  of 
the  algorithm  will  use  the  quotients  r/t  or  s/t,  assuming  that  t is  bigger  than 
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r and  s,  in  order  to  avoid  divisions  by  zero  or  by  near-zero  values.  If  r or  s is 
the  biggest,  it  is  recommended  to  adapt  the  algorithm  appropriately. 


VS  algorithm 

A '■=  Cd, 0,0 
for  i = 1 to  d 
B *=  Crf— i,t,0 
for  j = 1 to  i 

B :=  B ■ ( s/t ) + Cd—i'i—jj 
end  j 

A ■=  A ■ (r/t)  + B 

end  i 

p(r,  s,t)  = A ■ td. 

The  previous  algorithm  requires  ( d 2 + 5d)/2  multiplications  and  two  di- 
visions, and  so  is  significantly  faster  than  the  de  Casteljau  algorithm.  A 
backward  error  analysis  of  the  VS  algorithm  is  performed  in  the  following 
result. 

Theorem  4.  Consider  the  VS  algorithm  in  (15),  and  assume  that  4du  < 1. 
Then  the  computed  value  p(x)  = fl(p(x))  satisfies 

d i 

p(x)  (ie) 

i=0  j= 0 


where 


< 74d- 


(17) 


Proof:  The  VS  algorithm  consists  of  a Horner  type  algorithm  which  calculates 
and  a last  step  which  multiplies  by  td . From  (15)  we  can  write 


P{r,s,t) 

td 


EE/- 

t=0  j=o 


with  fd-i'i-jj  = . Since  by  (1)  fl(r/<)  = (r/f)(l  + 9\)  and  fl(s/<)  = 

(s/t)(  1 + 6*i ),  we  can  apply  Theorem  3.1  of  [7]  to  the  Horner  type  part  of  the 
VS-algorithm,  and  get 


fl 


where 


fd—i,i—j,j  — fd—i,i—j,j{l  T ^3d+l)- 


(18) 
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Finally,  in  the  last  step  we  have  to  multiply  by  td,  and  then  by  apply- 
ing d — 1 times  (1),  we  can  obtain  (16)  with  c.d-i,i-j,j  = = 

tdfd-i,i-j,j(l  + Od- 1).  Thus  by  (18)  and  (3) 


Cd-i,i-j,j  = + #3<i+l)(l  + &d- 1)  = tdfd-i,i-j,j(l  + 04  d) 

= d"  04 d)‘  1-1 

As  a consequence  of  Theorem  4,  and  applying  again  formula  (7),  we  can 
perform  a forward  error  analysis  of  the  VS  algorithm: 

Corollary  5.  Under  the  assumptions  of  Theorem  4,  the  computed  value 
p(r,s,t ) = fl(p(r,  s,  £))  of  the  VS  algorithm  satisfies 

|p(r,s,t)  — p(r,s,t)|  < ^4dCg(p(r,s,t)),  (19) 

where  B is  the  basis  used  in  the  modified  Bernstein-Bezier  representation. 

Although  in  Computer  Aided  Geometric  Design,  a bivariate  polynomial 
is  usually  stored  in  its  Bernstein  Bezier  form  (10)  (very  close  to  the  modified 
Bernstein-Bezier  representation  (15))  we  shall  compare  our  error  bounds  with 
those  obtained  by  evaluating  the  polynomial  in  Taylor  form  by  the  bivariate 
Horner  algorithm.  Given  the  triangle  T and  the  polynomial  p of  total  degree 
d defined  on  T,  let  u = x — xi,  v = y — yi,  where  (xi,yi)  are  the  cartesian 
coordinates  of  a point  of  T.  The  Taylor  form  of  p is  given  by 

d d—i 

p(u,v)  = J2Y^ai’iu'v3-  (20) 

i=0  j= o 


Bivariate  Horner  algorithm 


P :=  o-o, d 
for  i = 1 to  d 
A . — aitd — i 
for  j = 1 to  i 

A :=  A • u Qi—jtd—i 
end  j 
end  i 

p = A • v + A. 


We  observe  that  the  bivariate  Horner  algorithm  requires  (d2  + 3d) /2  mul- 
tiplications. 

The  following  result,  which  is  a consequence  of  Corollary  3.2  of  [7]  for  the 
backward  error  bound  and  of  formula  (7)  of  the  present  paper  for  the  forward 
error  bound,  provides  backward  and  forward  error  bounds  of  the  bivariate 
Horner  algorithm. 
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Theorem  6.  Consider  the  bivariate  Horner  algorithm  in  (20),  and  assume 
that  (2d  + 1)k  < 1.  Then  the  computed  value  p(u,v ) = fi(p(u,v))  satisfies 

d d—i 

p(u,  t>)  = °-i,ju'v3 , (21) 

i=0  j= o 


where 


and  satisfies 


| &i,j,k  ai,j,k\ 


< 72d+l, 


|p(m,^)  ~P(u,v) | < 72d+iCM(p(M,v)). 


(22) 

(23) 


§4.  Conclusions 

As  mentioned,  in  the  context  of  Computer  Aided  Geometric  Design,  polyno- 
mials are  usually  stored  in  the  Bernstein-Bezier  form  (10),  which  is  used  by 
the  de  Casteljau  algorithm.  Let  us  observe  that  the  coefficients  of  the  modified 
Bernstein-Bezier  form  (15)  used  by  the  VS  algorithm  are  related  with  those 
of (10)  by 


d\ 

cd-i,i-j,j  - (d  _ ^ 0 <j,i<d.  (24) 

The  conversion  from  the  Bernstein-Bezier  form  into  the  modified  Bernstein- 
Bezier  form  can  be  done  in  (d2  + 3d— 4)/2  multiplications.  In  [8]  the  algorithm 
composed  of  the  conversion  from  the  Bernstein-Bezier  form  into  the  modified 
Bernstein-Bezier  form  followed  by  the  VS  algorithm  was  called  the  VSC  al- 
gorithm. The  number  of  multiplications  and  divisions  required  by  the  VSC 
algorithm  is  d?  + id. 

We  have  seen  that  the  VS  algorithm  is  considerably  more  efficient  than  the 
de  Casteljau  algorithm.  On  the  other  hand,  the  bivariate  Horner  algorithm  is 
also  very  efficient  and  has  the  lowest  backward  error  bound,  as  one  can  deduce 
from  the  results  in  the  previous  section.  However  the  behaviour  with  respect 
to  the  forward  errors  is  different  as  we  shall  now  show. 

Given  the  bivariate  Bernstein  basis  B defined  on  a triangle  (see  (11)), 
since  the  barycentric  coordinates  satisfy  r + s -M  = 1 we  can  change  the 
variables  in  the  form  (r,  s,  t)  — (u,  v,  1 — u — v)  and  obtain  the  following 
expression  of  the  functions  in  B: 

Bfj(u,v)  — ^^ju'vj(l~u-v)d~'~3,  u,v  e [0,1], 

/ d\  _ d\ 

\hj)  ' (d  — i — 


where 
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With  the  same  change  of  variables,  the  functions  of  the  basis  B (used  in 
the  modified  Bernstein  form)  can  be  written  as 

Bf  j(u,v)  = u'yi (l  — u — v)d~l~\  u,v  € [0, 1], 

Thus,  these  functions  are  related  with  those  of  B by 


BijM  = 


We  see  that  the  basis  functions  of  B are  obtained  from  those  of  B by  a positive 
scaling.  Then,  using  the  condition  number  of  (5),  it  is  easy  to  prove  that 

CB(p(u,v))  = CB(p(u,v)),  Vp,  Vm,u  e [0,1].  (25) 


On  the  other  hand,  the  Taylor  form  uses  the  power  basis  M formed  by 
the  functions  uxv* , 0 < i < d,  0 < j < d — i,  u,  v € [0, 1].  By  formula  (100)  of 
[2],  the  functions  of  the  power  basis  M can  be  expressed  as 


= 


U V 


n~j n—k 


= 2.2- 


m 


( d ) 

=i  l=j  \i,j/ 


Since  the  coefficients  are  positive  and  the  basis  functions  are  nonnegative,  we 
can  deduce  from  Lemma  1 that  Cb(p{u,v))  < Cm{p(u,v)).  Therefore,  by 
(25),  we  conclude  for  every  polynomial  p(u,v), 


CB(p{u,v))  = CB(p{u,v))  <CM(p(u,v)),  u,v€[0,1].  (26) 


In  consequence,  although  the  bivariate  Horner  algorithm  provided  lower  back- 
ward error  bounds  than  de  Casteljau  and  VS  algorithms,  formula  (26)  shows 
that  these  algorithms  use  better  conditioned  bases,  and  this  fact  reduces  their 
corresponding  forward  error  bounds. 

In  conclusion,  the  VS  algorithm  has  more  advantages  than  de  Casteljau 
or  Horner  algorithms  in  this  context  due  to  the  following  properties.  First,  it 
uses  a basis  very  close  to  the  Bernstein  basis,  which  is  more  appropriate  in  the 
field  of  Computer  Aided  Geometric  Design  than  the  power  basis.  Secondly, 
the  bases  used  by  de  Casteljau  algorithm  and  the  VS  algorithm  are  also  better 
conditioned  than  the  power  basis  in  the  considered  domain,  this  last  property 
being  convenient  from  the  point  of  view  of  forward  error  analysis.  Finally,  the 
VS  algorithm  has  a higher  efficiency  than  the  de  Casteljau  algorithm. 
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